home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 4 / Apprentice-Release4.iso / Languages / RLaB 1.18c / testmatrix / trap2tri.r < prev    next >
Text File  |  1994-12-20  |  2KB  |  79 lines
  1. //-------------------------------------------------------------------//
  2.  
  3. // Synopsis:    Unitary reduction of trapezoidal matrix to triangular form.
  4.  
  5. // Syntax:    T = trap2tri ( L )
  6.  
  7. // Description:
  8.  
  9. //    Where L is an m-by-n lower trapezoidal matrix with m >= n,
  10. //    produces a unitary Q such that QL = [T; 0], where T is n-by-n
  11. //    and lower triangular. Q is a product of Householder
  12. //    transformations. 
  13.  
  14. //    Trap2tri returns a list with elements `q' and `t'.
  15.  
  16. //      Reference:
  17. //        G.H. Golub and C.F. Van Loan, Matrix Computations, second edition,
  18. //        Johns Hopkins University Press, Baltimore, Maryland, 1989.
  19. //        P5.2.5, p. 220.
  20.  
  21. //    This file is a translation of trap2tri.m from version 2.0 of
  22. //    "The Test Matrix Toolbox for Matlab", described in Numerical
  23. //    Analysis Report No. 237, December 1993, by N. J. Higham.
  24.  
  25. //-------------------------------------------------------------------//
  26.  
  27. trap2tri = function ( L )
  28. {
  29.   local (L)
  30.  
  31.   n = L.nr; r = L.nc;
  32.  
  33.   if (r > n || norm(L-tril(L),"1")) {
  34.     error("Matrix must be lower trapezoidal and m-by-n with m >= n.");
  35.   }
  36.  
  37.   Q = eye(n,n);        // To hold product of H.T.s
  38.  
  39.   if (r != n)
  40.   {
  41.     // Reduce nxr L =   r  [L1]  to lower triangular form: QL = [T].
  42.     //                 n-r [L2]                                 [0]
  43.  
  44.     for (j in r:1:-1)
  45.     {
  46.       // x is the vector to be reduced, which we overwrite with the H.T. vector.
  47.       x = L[j:n;j];
  48.       x[2:r-j+1] = zeros(r-j,1);    // These elts of column left unchanged.
  49.       s = norm(x,"2")*(sign(x[1]) + (x[1]==0));    // Modification for sign(1)=1.
  50.  
  51.       // Nothing to do if x is zero (or x=a*e_1, but we don't check for that).
  52.       if (s != 0)
  53.       {
  54.         x[1] = x[1] + s;
  55.         beta = s'*x[1];
  56.  
  57.         // Implicitly apply H.T. to pivot column.
  58.         // L(r+1:n,j) = zeros(n-r,1); // We throw these elts away at the end.
  59.         L[j;j] = -s;
  60.  
  61.         // Apply H.T. to rest of matrix.
  62.         if (j > 1)
  63.         {
  64.           y = x'*L[j:n; 1:j-1];
  65.           L[j:n; 1:j-1] = L[j:n; 1:j-1] - x*(y/beta);
  66.         }
  67.  
  68.         // Update H.T. product.
  69.         y = x'*Q[j:n;];
  70.         Q[j:n;] = Q[j:n;] - x*(y/beta);
  71.       }
  72.     }
  73.   }
  74.  
  75.   T = L[1:r;];    // Rows r+1:n have been zeroed out.
  76.  
  77.   return << q = Q; t = T >>;
  78. };
  79.